home *** CD-ROM | disk | FTP | other *** search
/ IRIX Base Documentation 1998 November / IRIX 6.5.2 Base Documentation November 1998.img / usr / share / catman / p_man / cat3 / complib / sggsvd.z / sggsvd
Text File  |  1998-10-30  |  9KB  |  265 lines
  1.  
  2.  
  3.  
  4. SSSSGGGGGGGGSSSSVVVVDDDD((((3333FFFF))))                                                          SSSSGGGGGGGGSSSSVVVVDDDD((((3333FFFF))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      SGGSVD - compute the generalized singular value decomposition (GSVD) of
  10.      an M-by-N real matrix A and P-by-N real matrix B
  11.  
  12. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  13.      SUBROUTINE SGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB,
  14.                         ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, IWORK, INFO
  15.                         )
  16.  
  17.          CHARACTER      JOBQ, JOBU, JOBV
  18.  
  19.          INTEGER        INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
  20.  
  21.          INTEGER        IWORK( * )
  22.  
  23.          REAL           A( LDA, * ), ALPHA( * ), B( LDB, * ), BETA( * ), Q(
  24.                         LDQ, * ), U( LDU, * ), V( LDV, * ), WORK( * )
  25.  
  26. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  27.      SGGSVD computes the generalized singular value decomposition (GSVD) of an
  28.      M-by-N real matrix A and P-by-N real matrix B:
  29.  
  30.          U'*A*Q = D1*( 0 R ),    V'*B*Q = D2*( 0 R )
  31.  
  32.      where U, V and Q are orthogonal matrices, and Z' is the transpose of Z.
  33.      Let K+L = the effective numerical rank of the matrix (A',B')', then R is
  34.      a K+L-by-K+L nonsingular upper triangular matrix, D1 and D2 are M-by-
  35.      (K+L) and P-by-(K+L) "diagonal" matrices and of the following structures,
  36.      respectively:
  37.  
  38.      If M-K-L >= 0,
  39.  
  40.                          K  L
  41.             D1 =     K ( I  0 )
  42.                      L ( 0  C )
  43.                  M-K-L ( 0  0 )
  44.  
  45.                        K  L
  46.             D2 =   L ( 0  S )
  47.                  P-L ( 0  0 )
  48.  
  49.                      N-K-L  K    L
  50.        ( 0 R ) = K (  0   R11  R12 )
  51.                  L (  0    0   R22 )
  52.  
  53.      where
  54.  
  55.        C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
  56.        S = diag( BETA(K+1),  ... , BETA(K+L) ),
  57.        C**2 + S**2 = I.
  58.  
  59.        R is stored in A(1:K+L,N-K-L+1:N) on exit.
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. SSSSGGGGGGGGSSSSVVVVDDDD((((3333FFFF))))                                                          SSSSGGGGGGGGSSSSVVVVDDDD((((3333FFFF))))
  71.  
  72.  
  73.  
  74.      If M-K-L < 0,
  75.  
  76.                        K M-K K+L-M
  77.             D1 =   K ( I  0    0   )
  78.                  M-K ( 0  C    0   )
  79.  
  80.                          K M-K K+L-M
  81.             D2 =   M-K ( 0  S    0  )
  82.                  K+L-M ( 0  0    I  )
  83.                    P-L ( 0  0    0  )
  84.  
  85.                         N-K-L  K   M-K  K+L-M
  86.        ( 0 R ) =     K ( 0    R11  R12  R13  )
  87.                    M-K ( 0     0   R22  R23  )
  88.                  K+L-M ( 0     0    0   R33  )
  89.  
  90.      where
  91.  
  92.        C = diag( ALPHA(K+1), ... , ALPHA(M) ),
  93.        S = diag( BETA(K+1),  ... , BETA(M) ),
  94.        C**2 + S**2 = I.
  95.  
  96.        (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
  97.        ( 0  R22 R23 )
  98.        in B(M-K+1:L,N+M-K-L+1:N) on exit.
  99.  
  100.      The routine computes C, S, R, and optionally the orthogonal
  101.      transformation matrices U, V and Q.
  102.  
  103.      In particular, if B is an N-by-N nonsingular matrix, then the GSVD of A
  104.      and B implicitly gives the SVD of A*inv(B):
  105.                           A*inv(B) = U*(D1*inv(D2))*V'.
  106.      If ( A',B')' has orthonormal columns, then the GSVD of A and B is also
  107.      equal to the CS decomposition of A and B. Furthermore, the GSVD can be
  108.      used to derive the solution of the eigenvalue problem:
  109.                           A'*A x = lambda* B'*B x.
  110.      In some literature, the GSVD of A and B is presented in the form
  111.                       U'*A*X = ( 0 D1 ),   V'*B*X = ( 0 D2 )
  112.      where U and V are orthogonal and X is nonsingular, D1 and D2 are
  113.      ``diagonal''.  The former GSVD form can be converted to the latter form
  114.      by taking the nonsingular matrix X as
  115.  
  116.                           X = Q*( I   0    )
  117.                                 ( 0 inv(R) ).
  118.  
  119.  
  120. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  121.      JOBU    (input) CHARACTER*1
  122.              = 'U':  Orthogonal matrix U is computed;
  123.              = 'N':  U is not computed.
  124.  
  125.  
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. SSSSGGGGGGGGSSSSVVVVDDDD((((3333FFFF))))                                                          SSSSGGGGGGGGSSSSVVVVDDDD((((3333FFFF))))
  137.  
  138.  
  139.  
  140.      JOBV    (input) CHARACTER*1
  141.              = 'V':  Orthogonal matrix V is computed;
  142.              = 'N':  V is not computed.
  143.  
  144.      JOBQ    (input) CHARACTER*1
  145.              = 'Q':  Orthogonal matrix Q is computed;
  146.              = 'N':  Q is not computed.
  147.  
  148.      M       (input) INTEGER
  149.              The number of rows of the matrix A.  M >= 0.
  150.  
  151.      N       (input) INTEGER
  152.              The number of columns of the matrices A and B.  N >= 0.
  153.  
  154.      P       (input) INTEGER
  155.              The number of rows of the matrix B.  P >= 0.
  156.  
  157.      K       (output) INTEGER
  158.              L       (output) INTEGER On exit, K and L specify the dimension
  159.              of the subblocks described in the Purpose section.  K + L =
  160.              effective numerical rank of (A',B')'.
  161.  
  162.      A       (input/output) REAL array, dimension (LDA,N)
  163.              On entry, the M-by-N matrix A.  On exit, A contains the
  164.              triangular matrix R, or part of R.  See Purpose for details.
  165.  
  166.      LDA     (input) INTEGER
  167.              The leading dimension of the array A. LDA >= max(1,M).
  168.  
  169.      B       (input/output) REAL array, dimension (LDB,N)
  170.              On entry, the P-by-N matrix B.  On exit, B contains the
  171.              triangular matrix R if M-K-L < 0.  See Purpose for details.
  172.  
  173.      LDB     (input) INTEGER
  174.              The leading dimension of the array B. LDA >= max(1,P).
  175.  
  176.      ALPHA   (output) REAL array, dimension (N)
  177.              BETA    (output) REAL array, dimension (N) On exit, ALPHA and
  178.              BETA contain the generalized singular value pairs of A and B;
  179.              ALPHA(1:K) = 1,
  180.              BETA(1:K)  = 0, and if M-K-L >= 0, ALPHA(K+1:K+L) = C,
  181.              BETA(K+1:K+L)  = S, or if M-K-L < 0, ALPHA(K+1:M)=C,
  182.              ALPHA(M+1:K+L)=0
  183.              BETA(K+1:M) =S, BETA(M+1:K+L) =1 and ALPHA(K+L+1:N) = 0
  184.              BETA(K+L+1:N)  = 0
  185.  
  186.      U       (output) REAL array, dimension (LDU,M)
  187.              If JOBU = 'U', U contains the M-by-M orthogonal matrix U.  If
  188.              JOBU = 'N', U is not referenced.
  189.  
  190.  
  191.  
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. SSSSGGGGGGGGSSSSVVVVDDDD((((3333FFFF))))                                                          SSSSGGGGGGGGSSSSVVVVDDDD((((3333FFFF))))
  203.  
  204.  
  205.  
  206.      LDU     (input) INTEGER
  207.              The leading dimension of the array U. LDU >= max(1,M) if JOBU =
  208.              'U'; LDU >= 1 otherwise.
  209.  
  210.      V       (output) REAL array, dimension (LDV,P)
  211.              If JOBV = 'V', V contains the P-by-P orthogonal matrix V.  If
  212.              JOBV = 'N', V is not referenced.
  213.  
  214.      LDV     (input) INTEGER
  215.              The leading dimension of the array V. LDV >= max(1,P) if JOBV =
  216.              'V'; LDV >= 1 otherwise.
  217.  
  218.      Q       (output) REAL array, dimension (LDQ,N)
  219.              If JOBQ = 'Q', Q contains the N-by-N orthogonal matrix Q.  If
  220.              JOBQ = 'N', Q is not referenced.
  221.  
  222.      LDQ     (input) INTEGER
  223.              The leading dimension of the array Q. LDQ >= max(1,N) if JOBQ =
  224.              'Q'; LDQ >= 1 otherwise.
  225.  
  226.      WORK    (workspace) REAL array,
  227.              dimension (max(3*N,M,P)+N)
  228.  
  229.      IWORK   (workspace) INTEGER array, dimension (N)
  230.  
  231.      INFO    (output)INTEGER
  232.              = 0:  successful exit
  233.              < 0:  if INFO = -i, the i-th argument had an illegal value.
  234.              > 0:  if INFO = 1, the Jacobi-type procedure failed to converge.
  235.              For further details, see subroutine STGSJA.
  236.  
  237. PPPPAAAARRRRAAAAMMMMEEEETTTTEEEERRRRSSSS
  238.      TOLA    REAL
  239.              TOLB    REAL TOLA and TOLB are the thresholds to determine the
  240.              effective rank of (A',B')'. Generally, they are set to TOLA =
  241.              MAX(M,N)*norm(A)*MACHEPS, TOLB = MAX(P,N)*norm(B)*MACHEPS.  The
  242.              size of TOLA and TOLB may affect the size of backward errors of
  243.              the decomposition.
  244.  
  245.  
  246.  
  247.  
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.